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Abstract. We consider an algorithm called FEMWARP for warping tetrahedral finite element meshes that 
computes the warping using the finite element method itself. The algorithm takes as input a two- or three-dimensional 
domain defined by a boundary mesh (segments in one dimension or triangles in two dimensions) that has a volume 
mesh (triangles in two dimensions or tetrahedra in three dimensions) in its interior. It also takes as input a prescribed 
movement of the boundary mesh. It computes as output updated positions of the nodes of the volume mesh. The 
first step of the algorithm is to determine from the initial mesh a set of local weights for each interior node that 
describes each interior node in terms of the positions of its neighbors. These weights are computed using a finite 
element stiffness matrix. After a boundary transformation is applied, a linear system of equations based upon the 
weights is solved to determine the final positions of the interior nodes. 

Our main concern is when this algorithm reverses elements. We analyze the causes for this undesirable behavior 
and propose techniques to make the method more robust against reversals. Among the methods includes combining 
FEMWARP with an optimization-based untangler. 

Finally, we use FEMWARP to study the movement of the canine heart. 

Key words, moving meshes, adaptation, finite element method, optimization-based mesh untangling, tetrahedral 
meshes, unstructured mesh generation, cardiology 
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1. Introduction. Moving meshes arise in cardiology, computer graphics, animation, and crash 
simulation, among other applications in science and engineering. With moving meshes, the mesh 
is updated at each step in time due to a moving domain boundary, thus resulting in potentially 
drastically varying mesh quality from step to step. One problem that can occur at each timestep 
is element reversal. "Reversal" means that the element changes orientation. In two dimensions, 
this means that the nodes of an element are clockwise when they ought to be counterclockwise, 
and in three dimensions it means that they violate the right-hand rule. We focus on mesh-warping 
procedures that avoid reversals as much as possible. 

It is well-known that poor quality elements affect the stability, convergence, and accuracy of 
finite element and other solvers because they result in poorly-conditioned stiffness matrices and 
poor solution approximation j^T] . If well-shaped elements are not the result of updating the mesh 
boundary, the mesh quality must be improved by topological or geometrical means after each time 
step. 

Research has shown that mesh smoothing (or r-refinement) methods can be applied to improve 
the quality of a mesh. These methods adjust the positions of the vertices in the mesh while preserving 
its topology. 

Laplacian smoothing is the most popular method for node-based mesh smoothing. In an iterative 
manner, it repositions the vertices of the mesh by moving each interior node to the geometric center 
of its neighbors. It is often used because it is computationally inexpensive and is very easy to 
implement. However, the method has several undesirable properties. One of them is that the 
method is not guaranteed to work, i.e., sometimes it reverses mesh elements. A second drawback is 
that if the algorithm is not run to convergence, the resulting mesh depends upon the order in which 
the nodes are smoothed. 
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A related type of smoothing, namely Winslow smoothing, is more resistant to mesh folding due 
to the requirement that the logical variables be harmonic functions. This method was introduced for 
structured meshes by Winslow in '32' . In , Knupp extends this idea to unstructured quadrilateral 
meshes; he used a finite-difference method for his derivative computations. Others have used the 
finite element method for the derivative computations on unstructured meshes; for example, Knupp 
cites a proprietary DOE paper by Tipton to which the authors do not have access. (See the references 
at the end of 

Many others have created extensions to Winslow's method. One type of extension has been 
to generate the meshes using a variational approach; for example, see Brackbill and Saltzman fJli 
Russell and co-workers |B] and [21 , and Thompson et al. . Dvinsky |Sj uses the theory of harmonic 
functions to generate a harmonic map between the physical domain and the logical domain in order 
to generate the mesh. Meshes that have been generated according to this method possess regularity 
and smoothness properties. Li et al. 22 also introduced a moving mesh method based on harmonic 
maps; their method has two parts: a solution algorithm and a mesh-redistribution algorithm. In 
doing so, they are able to obtain the desirable properties of both the h-method and r-mcthod for 
finite elements. 

Other, more accurate methods for r-refinement are possible. Most of these methods are based 
upon optimization. Optimization-based methods are used with the goal of guaranteeing an improve- 
ment in the mesh quality by minimizing a particular mesh quality metric. Their main drawback, 
however, is their computational expense. Examples of optimization-based methods for r-refinement 
can be found in the following papers: [El, E], [El, 0, [H, IB, [2], ESI, 0, ED, and ES]. 
For a theory of algebraic mesh quality metrics see EO] • 

To address the above issues. Baker El developed a three-step method for the metamorphosis 
of tetrahedral meshes. Each cycle of the method involves a combination of r-refinement, mesh 
coarsening, and mesh enrichment to adapt the mesh. The first step in the cycle is to move the 
interior nodes as far as possible using r-refinement while avoiding element reversal. The second step 
is to remove the poorly-shaped elements in the mesh using mesh coarsening. The final step is the 
addition of elements to improve the mesh quality by mesh refinement. This dynamic procedure was 
shown to be more cost-effective than just r-refinement. One disadvantage of this technique is that it 
comes with no theoretical guarantees as it is difficult to analyze because each cycle is a combination 
of three very different techniques. 

We study a different mesh warping problem where the connectivity of the mesh is not allowed to 
change, which is important for some applications. The problem that we address is as follows: Given 
a three-dimensional domain, bounded by a triangulated surface mesh, and given an interior volume 
mesh composed of unstructured tetrahedra, suppose the triangulated surface mesh is displaced. Is 
there an algorithm to move the nodes of the volume mesh so that it continues to conform to the 
surface mesh and to be a good quality mesh? In addition, we study the analogous problem in 
two dimensions. Seeking simplicity, preservation of the mesh's combinatorial structure, theoretical 
insight, and low computational expense, we propose a technique called FEMWARP and several 
variants to be described below. 

The FEMWARP algorithm has already been considered in the previous literature, e.g., by 
Baker EI, although he rejected it in favor of a method based on linear elasticity. It is shown, 
however, in that there appear to be few advantages of linear elasticity over FEMWARP; whereas, 
a significant disadvantage is that the linear elasticity matrix problem is three times larger (in 3D) 
than FEMWARP 's (although FEMWARP must solve the smaller linear system three times). 

FEMWARP is described in Section [2 and a generalization of it called the "linear weighted 
Laplacian smoothing" (LWLS) framework is considered in Section|3| There are two main advantages 
to methods within the LWLS framework. First, if a continuous deformation of the boundary is given, 
then methods within the framework are valid for computing the resulting trajectory that specifies 
the movement of the interior nodes. In addition, these trajectories will be continuous. This is vital 
for some applications where continuity of motion is required. A second big advantage is that sparse 
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matrix algorithms may be used to solve (|2.2|l . The sparsity structure is apparent, since, on average, 
an interior node has six neighbors in 2D, whereas a typical 2D mesh may have thousands of nodes. In 
Section^ we explain the relationship between traditional Laplacian smoothing and LWLS methods. 

The principal failure mode for FEMWARP is that it can reverse elements. The causes of 
element reversal are covered in Section |S1 Techniques to prevent some reversals, including small- 
step FEMWARP and mesh refinement, are also covered in that section. Another technique to avoid 
reversals based on the Opt-MS mesh untangler is covered in Sectional In Section [7| we test our 
algorithms on several types of deformations of three-dimensional meshes. In Section |S1 we apply 
FEMWARP to study the motion of the beating canine heart under normal conditions. 

2. The FEMWARP algorithm. The first step of the FEMWARP algorithm is to express 
the coordinates of each interior node of the initial mesh as a linear combination of its neighbors. Let 
a triangular or tetrahedral mesh be given for the domain Q in two or three dimensions. Let b and 
m represent the numbers of boundary and interior nodes, respectively. Form the (m -I- 6) x (m -I- b) 
stiffness matrix A based on piecewise linear finite elements defined on the initial mesh for the 
boundary value problem 

Au — on fl. 

It is well-known |18| that this matrix is determined as follows. Let (f>i be the continuous piecewise 
linear function (where the pieces of linearity are given by the triangulation) such that 4>i{xi) = 1, 
where Xi is the i*'' node of the mesh, and 4>i{xj) = 0, where Xj is any other node in the mesh (j 7^ i). 
Define for each j, j — 1, . . . , m -I- 6 

Jn 

This matrix will be sparse and symmetric positive semidefinite. Its nonzero entries correspond to 
pairs of neighboring nodes in the mesh. 

Next, let Ai denote the m x m submatrix of A whose rows and columns are indexed by interior 
nodes, and let Ab denote the m x b submatrix of A whose rows are indexed by interior nodes 
and whose columns are indexed by boundary nodes. Let x be the (rn + b) vector consisting of x- 
coordinates of the nodes of the initial mesh, where we assume that the interior nodes are numbered 
first. Then it follows from well-known theory that [Aj, Agjx = because any linear function of the 
coordinates is in the null-space of the discretized Laplacian operator. For the same reason, a similar 
identity holds for the y- and ^-coordinates. An equivalent way to write this equation is 

Aixi^-AbXb- (2.1) 

If we divide each row of [Aj, Ab] by the diagonal element in that row, we obtain a linear system 
whose diagonal entries are I's and whose row sums are O's. This means that the [A/,j4b], thus 
scaled, expresses each interior nodal coordinate as an affine combination of the neighboring nodal 
coordinates. 

The formation of Ai and is the first step of the FEMWARP method. Consider now the 
application of a user-supplied transformation to the boundary of the mesh. We denote the new 
positions of the boundary nodes by [xb, Vb] in two dimensions or [x_b, ys, zb] in three. 

The final step is to solve a linear system of equations similar to 12.1|l for the new coordinates of 
the interior nodes. We solve H2.2(l for [xi,yi] : 

Ai[xi,yi]^ -AB[xB,yB]- (2.2) 
or the analog in three dimensions. This concludes the description of FEMWARP. 
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3. The LWLS framework. The FEMWARP method can be generahzed as follows. One 
produces, using some method, a family of weights Wij for each ordered pair of nodes (i, j) such that 
j is a neighbor of i and such that i is interior. Denote by N{i) the set of neighbors of i, and let 
the coordinates of node i be (xi,yi) in 2D or {xi,yi,Zi) in 3D. We want these weights to express 
interior nodes as affine combinations of their neighbors, which imposes the following constraints on 
the Wij's: For each i indexing an interior node, 

J2 = 1, 

j<£N{2) 

and in three dimensions, there is a fourth equation for z-coordinates. Let [A/jvls] be to x (m + b) 
matrix with I's on the diagonal and —Wij in position {i,j) whenever j £ The above three 

equations may be rewritten in terms of Aj and as follows: 

Aid + AbCb = 0, 
Aixi + Abxb = 0, 
Ajyi + AbVb = 0, 

where e/ and cb are vectors of all I's whose length is equal to the number of interior and boundary 
nodes respectively, xj and xb are the vectors of x-coordinates of mesh nodes in the interior and 
boundary respectively, and similar for yB and yj. 

Once these weights are derived, then, as in FEMWARP, the user-supplied deformation may be 
applied to the boundary nodes yielding deformed boundary coordinates (xb, ys) in two dimensions. 
Finally, the new position of the interior nodes is computed via 

Ai[xi,yi] = -AB[xB,yB]- 

We will say that any method following this framework is a linearly weighted Laplacian smooth- 
ing (LWLS) method. The precise relationship between these methods and traditional Laplacian 
smoothing is the subject of Sectional The FEMWARP method was not written exactly this way 
since the diagonal elements in [A/,^!^] for FEMWARP are not rescaled to be equal to 1. Usually 
for FEMWARP it is preferable to preserve the symmetry of Aj and therefore omit the scaling step. 
This omission has no effect on the final answer. 

One concern with the weights in FEMWARP is that, in general, they are not always nonnegative. 
Sometimes it is desirable to have all nonnegative weights. Geometrically, this means that the interior 
nodes are expressed as convex (rather than merely affine) combinations of neighbors. Nonnegative 
weights are useful if the weight matrix is used to interpolate properties besides nodal coordinates, 
and it is mandatory that no extrapolation occur (e.g., because a nodal value that strays outside the 
range of boundary data would violate some kind of physical constraint). 

In |28| . the first author experimented with other weighting schemes including a method called 
LBWARP in which the weights Wij are obtained by solving a log barrier optimization problem for 
each i: 

max EiGJVwlog^'y 
subject to J2jeN(t) '^ij = 1 
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The rationale for this optimization formulation is to ensure that each Wij is positive, and, more 
strongly, that it is bounded well away from (since the logarithmic term tends to —oo if any of the 
Wij 's approaches zero) . This problem is convex, and there is an efficient algorithm to find the unique 
global minimum. It can be proved 28 that the w^ 's that solve this optimization problem have a 
positive lower bound that depends only on the worst-case aspect ratio of triangles or tetrahedra in 
the mesh. 

A further important property that a method in the LWLS family ought to satisfy is that Aj 
should be nonsingular to ensure a unique solution for the new coordinates. For FEMWARP, the 
nonsingularity of Aj follows from well-known theory of finite element analysis. For LBWARP, 
the nonsingularity follows because Aj is an irreducible and diagonally dominant matrix with strict 
inequality for at least one row. 

Because the mesh is connected and a positive weight is associated with each edge, we see that 
Aj is irreducible. Also, Aj is digonally dominant, as its diagonal entries are 1, and its off-diagonal 
entries are nonpositive and sum to a number in [—1,0] in each row. Since there is at least one 
interior node adjacent to a boundary node, Ajei ^ 0, where e/ is a vector of all I's of length m. 
Therefore, Ai satisfies the definition of an irreducibly diagonally dominant matrix by the above 
argument. Thus, Aj is invertible PHI Theorem 1.8]. 

Experiments in indicate that weights obtained from LBWARP generally do not perform as 
well as FEMWARP for the mesh warping application in terms of resisting element reversal. 

One useful property that holds for any method in the LWLS framework including FEMWARP 
and LBWARP is that the method is exact for affine transformations. Let us state this as a theorem. 
The theorem is stated for the two-dimensional case, and it extends in the obvious way to three 
dimensions. 

Theorem 3.1. Let As and Aj be generated using a method from the LWLS framework, and 
assume Aj is nonsingular. Let [xbtVb] be the user- specified deformed coordinates of the boundary. 
Suppose there exists a 2 x 2 nonsingular matrix L and 2-vector v such that for each j Cz B, 




Let [xi , yi] be the deformed interior coordinates computed by the method. Then for each i Cz I , 




Proof. The positions of the interior nodes in the deformed mesh are given by 

[xi,yi] = -Aj^AB{[xB,yB]L^ + eBV^) (3.1) 

where, as above, xi,yi are column vectors composed of the x- and y-coordinates of the interior 
nodes respectively and xb, yB are the corresponding vectors for boundary nodes, and finally eg is 
vector of all I's of length b. 

In order to show that affine mappings yield exact results with any algorithm within the frame- 
work, we want to show that 13. 1|) is the same as: 

[xi,yi] = [xi,yi]L'^ + eiv^. (3.2) 

Observe that the equivalence of H3.1|) and H3.2I) would follow immediately from: 

Ai{[xi,yi]L'^ + eiv'^) = -^^([xs, ys]^^ + esw^). (3.3) 

Thus, it remains to check that (|3.3I) holds. 
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Because the weights for each interior node sum to 1, Ajcj + Ageg = as noted above. Hence 
{AjCj + ABeB)v'^ — 0. Also, because [xj, yj] and [xg, yg] denote the original positions of the nodes, 
we know that Ai[xi,yi] + ^^[xs, ys] = 0. So {Ai[xi,yi] + AB[xB,yB])L'^ = 0. 

Putting these together, we see that 

{Ai[xi, yi] + AB[xB,yB])L'^ + {Aiei + ABei)v^ = 0. (3.4) 
Therefore, holds, and the lemma is proven. □ 

4. Relationship to Laplacian smoothing. In this section we explain the fairly direct con- 
nection between Laplacian smoothing and algorithms from the LWLS framework. We start with the 
following theorem, which is a consequence of well-known results in the previous literature. 

Theorem 4.1. Gauss-Seidel iteration will converge for solving the linear system 

Ai[xi,yi] ^ -AB[xB,yB] 

that is produced by either the LBWARP or FEMWARP algorithm. 

Proof. In the case of FEMWARP, Ai is symmetric and positive definite hence Gauss-Seidel is 
convergent Theorem 10.1.2]. In the case of LBWARP, Ai is irreducibly diagonally dominant 
hence Gauss-Seidel is convergent [HOI Theorem 3.4]. □ 

The Gauss-Seidel algorithm applied to these matrices corresponds to iteratively replacing the 
coordinates of each interior vertex with a weighted average of the coordinates of its neighbors. The 
weights come from the entries of Ai and Ab- In the case of FEMWARP, some of the weights may 
be negative. 

Ordinary laplacian smoothing, as mentioned in Section^ corresponds to the iterative process of 
replacing each mesh node with the centroid of its neighbors. Thus, it can be regarded as a member 
of the LWLS framework in which the weight of each neighbor of node i is l/|A^(i)|, where N{i) is 
the set of neighbors of node i. 

5. Element reversal and small-step FEMWARP. Sometimes FEMWARP fails to yield a 
valid triangulation because it reverses elements. The purpose of this section is to explore the causes 
for reversal and propose some workarounds. The discussion of workarounds continues into the next 
section. 

Let be the original polygonal or polyhedral domain, and let be the domain whose boundary 
is given by the user-specified deformation of the boundary nodes of fi, i.e., by the nodes at coor- 
dinates {xB,yB)- Assume that this user-specified deformation is not self- intersecting and preserves 
orientation. 

Let (p be the mapping from O to computed by FEMWARP. In other words, interpolate the 
interior nodal deformations linearly over the elements to arrive at a continuous function on the whole 
domain. (In the case that FEMWARP fails, parts of 0(ri) may protrude outside of f2.) 

Let (jj* be the mapping that is obtained from the exact (continuum) Laplacian. In other words, 
solve the boundary value problems 

Ax = on n, 
Ay = on Q, 

X — xb on dQ, 

y ^ Vb on dfl 

and define (/)*(a;, y) = {x{x,y),y(x,y)). Let us call this warping algorithm "continuum FEMWARP" . 

Finally, let 0+ be the mapping that is obtained by linear interpolation over the elements of (f)* 
evaluated at nodes. Thus, 0^ is intermediate between (j) and (f>* in the sense that for (f)~^ we use the 
exact solution to the continuum problem only at nodal points and use interpolation elsewhere. 

There are three possible reasons that FEMWARP could fail: 
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Fig. 5.1. All meshes in this section were generated by the Triangle mesh generator; this is an example of a mesh 
used herein. 



1. Mapping 0* might have reversals. 

2. Mapping 0+ might have reversals even though 0* has none. 

3. Mapping might have reversals even though (j)^ has none. 

Let us call these Type 1, Type 2, and Type 3 failures. For Type 1 failure, "reversal" means the 
existence of a point x £ SI such that V4>{x) has a nonpositive determinant. For the second and third 
types, "reversal" means that a triangle is reversed. Type 1 reversals are caused by the boundary 
deformation alone and are not related to the mesh. Type 2 reversals are due to continuous versus 
discrete representation of 0* , and Type 3 reversals can be analyzed using traditional error estimates 
for the finite element method. 

Let us start with Type 1 reversals. It is difhcult to characterize inputs for which (jj* will have 
reversals. In j28| . a sufficient condition is given for two-dimensional domains to ensure that 0* will 
be an invertible function, but the condition is unrealistically stringent and is nontrivial to check in 
practice. 

Rather than presenting the theorem from |2H|) we choose to present a series of examples and 
discussion on how to avoid reversals. The geometry for most of the examples in this section is 
a two-dimensional annulus with outer radius 1 and inner radius r < 1. Meshes used for tests in 
this section were generated by Triangle, a two-dimensional quality mesh generation package j26| . 
A typical mesh that occurs in some of the experiments is depicted in Fig. 15.11 this mesh has 1238 
elements, 694 nodes, and a maximum side length of 0.1137. Its inner radius is 0.5. 

The boundary transformations applied in this section usually consist of two motions: First, 
the inner circular boundary is moved radially outward ( "compression" ) to new radius s such that 
r < s < 1. Second, a rotation of magnitude 9 is applied to the outer boundary. The relevant Laplace 
equations determining (f>* can be solved in closed form to yield 

(f)* (x, y) = {Ax + By,-Bx + Ay) 

where A = a + h/{x^ + y^) and B = c + djix^ -I- y^), and a,b,c,d are constants determined by 
the boundary conditions. In particular, to match the boundary conditions just described, one must 
satisfy the equations a + b ~ cos 6, c + d ^ sin 6, a + b/r'^ = s/r, c + d/r^ = 0. These equations are 
uniquely solved by choosing 

r , n {cos(0) — rs, rs — cos(6'), sin(0), — sin(0)} 

{a,b,c,d} — . 

1 — 

This function (j)* is invertible, i.e., avoid reversals, provided the determinant of its Jacobian is 
always positive. This determinant may be computed in closed form: det(V(/)*) = 0^ + 0^ — (6^ -I- 
d^)/{x'^ + y'^Y ■ This quantity is minimized when x"^ ^y"^ — r^. Therefore, reversals of Type 1 occur 
if and only if r'^{a? + c^) < 6^ + d^ . Substituting the above formulas for a,b,c,d yields the result 
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that reversals occur if and only if 

2rcos(6') - r'^s - s < 0. 

For example, if r = s = 0.5 (no compression), then reversals occur when cos(^) < .625, i.e., 
9 > 51.4...°. If r = .5 while s = .75, then reversals occur when 9 > 20.4...°. 

We tested the FEMWARP algorithm on the cases described above with a mesh for the annuluar 
region as discussed earlier. We used a mesh with inner radius r = 0.5. This particular mesh 
contained 10,950 triangles with maximum side length of 0.039. For s = 0.5, when we selected 
9 = 51°, FEMWARP ran on this mesh without reversals, whereas 6 — 52° caused reversals. As 
mentioned in the previous paragraph, 9 w 51.4° is the cutoff for Type 1 reversals. When r = 0.5, 
s = 0.75, FEMWARP succeeded for 9 = 22° but failed for 9 = 23°, again, very close to the cutoff 
for Type 1 reversals. 

The point of the experiments in the last paragraph is that for a reasonably refined and reason- 
ably high-quality mesh, most FEMWARP reversals seem to be Type 1 reversals. In other words, 
FEMWARP fails when continuum FEMWARP fails or is close to failure. Other experiments not 
reported here seem to confirm this point. Therefore, in order to extend the range of deformations 
that can be handled by FEMWARP, the best strategy is to come up with a way to avoid Type 1 
reversals. 

One simple method to avoid Type 1 reversals is to take several smaller steps instead of one big 
step. For example, suppose {x'g,y'g) are positions for the boundary nodes intermediate between 
their initial positions and their final positions {xB,yB)- Then one could define a two-step continuum 
FEMWARP as follows. Solve 

Ax' =0 on fi. 

Ay' =0 on n, 
x' = x'g on dQ, 

y' = tj'b on an 

for x' and y' to determine a mapping : ^ Q! given by {x, y) ^ {x' , y') (where f2' is the domain 
bounded by (^3,^3)) followed by 

Ax = on fi'. 
Ay = on Q', 

X = xb on dO.' , 

y = yB on dCl' 

fov X, y to obtain a map ^2- Finally, (f)* = <p2°(f>i- 

The idea in the previous paragraph can be extended to more steps with smaller increments. 
The limiting case of an infinite number of infinitesimal steps yields an algorithm that we call 
"infinitesimal-step continuum FEMWARP" . This algorithm corresponds to solving a time-dependent 
system of partial differential equations, but it is difficult even to write down the system that de- 
scribes this limit because it requires notation for inverting systems of Laplacians. We suspect that 
infinitesimal-step continuum FEMWARP will never suffer from reversals as long as the boundary 
does not get tangled. 

In the case of the annulus, it is possible again to write down infinitesimal-step continuum 
FEMWARP in closed form. For simplicity, let us assume r = s so that the only deformation is 
the rotation of the outer boundary. Assume this rotation is broken up into infinitesimally small 
rotations. (Another choice would be to connect the initial positions to the final positions with line 
segments, and break up the boundary motion as infinitesimal increments along the line segments. 
This way to obtain a continuous boundary motion is undesirable, however, because for a sufficiently 
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large rotation, the line segments would cut through the inner boundary of the annulus and hence 
cause tangling of the boundaries.) 

With the setup described in the last paragraph, the deformation for an outer rotation of 9 
computed by infinitesimal-step continuum FEMWARP maps a point at initial position p(cos 0, sin(f)) 
(r < p < 1) to p(cos((?!)+a), sin((/)-|-Q!)), where a = (1 — r^/p^)0/(l — r^). This map is clearly bijective 
for any value of 9; it corresponds to rotating each concentric circle of the annulus by an amount that 
interpolates between (when p = r) and 9 (when p = 1). 

Thus, by using small-step FEMWARP with sufficiently small steps, we can essentially eliminate 
Type 1 failures. Small-step FEMWARP preserves the attractive property of FEMWARP that it is 
exact for affine maps, as long as all the intermediate steps are also affine. Unfortunately, it loses 
the attractive property that only one coefficient matrix for solving the linear system needs to be 
factored. Small-step FEMWARP requires the solution of a different coefficient matrix for each step. 
This drawback is partly ameliorated by the fact that even though the matrices are different, they 
have the same nonzero pattern, and hence the symbolic phase of sparse direct solution may be 
reused. If instead an iterative method is being used to solve the mesh warping equations, then the 
sparsity pattern may be reused in the preconditioner. In addition, the factored coefficient matrix at 
step tk can be used as a preconditioner for an iterative method at step tk+i- 

Elimination of Type 1 failures means that the mapping function <j)* has no reversals in the sense 
that the determinant of its Jacobian is positive everywhere; equivalently, it does not reverse any 
infinitesimally small triangles. A Type 2 failure occurs because the triangles in the mesh have finite 
(non- infinitesimal) size and hence can still be reversed by 0"*". The following theorem characterizes 
when this can happen. 

Theorem 5.1. Suppose that f : fl —> is bijective, orientation-preserving and on Q, with 
V/ nonsingular. Let T he a triangle in the mesh with vertices {vi,V2,v^} , and let T' be the triangle 
whose vertices are {/(I'l), /(W2), /(t'a)}- If below holds, then T' is not reversed. 

Proof. Recall that triangle T with vertices {wi,W2,W3} is positively oriented if and only if 
det{A) > 0, where 

A= {V2- Vi,V3 - Vi). 

In order to analyze the analogous quantity for {/(wi), /(i'2), /(^s)}, we start with the following 
algebra, which invokes the fundamental theorem of calculus twice: 

f{v2)-fivi)^ f Wfi{l-t)v,+tV2){v2-Vi)dt 
JO 

1 



Vf{{l^t)vi+tV2)dt] {V2-Vi) 







V/(l'l) + / [V/((l - t)v, + tV2) - V/(l>i)] dt {V2 - I'l) 







V/(i;i) + 







V^t/((1 - S)VX + S((l - t)vx + tV2)){v2 ~ vi)ds 







dt {V2 - Vl) 



V/(wi)(u2 - Vl) + ei 



where ||ei|| < Mh'^, where h is the maximum side length of T (an upper bound on \\v2 — vi\\) and 
M is an upper bound on ||V^/|| in the triangle. Similarly, 

/(«3) - f{vi) = V/(t;i)(-y3 - Vl) + 62, 

where again |je2|| < Mh'^. Therefore, 

{f{v2) - f{vi), f{v3) ~ f{vi)) - V/(«i)A + E (5.2) 
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where A is as above and ||i?||2 < \/2Mh'^. Observe that \7f{vi)A has positive determinant by as- 
sumption. Therefore, the left-hand side can have negative determinant only if i? is a sufficiently large 
perturbation to change the determinant sign. If E is such a large perturbation, then by the continuity 
of the determinant, there is a perturbation E' no larger than E such that the det{W f (vi) A + E') = 0, 
i.e., V/(ui)A + £" is singular. Furthermore, < \/2AIh^. By Theorem 2.5.3 of ^Hli this means 
that 

V2Mh^ > C7„in(V/(«i)yl) > C7„in(V/(«i))a„,in(A), (5.3) 

where (Jmin(^) and <7rna.x{A) denote the smallest and largest singular values of A, respectively. It 
follows from the equation AA~^ = / that that the columns of A~-^ are parallel to the altitude 
segments of triangle T perpendicular to viv^ and viV2 respectively, but scaled so that their lengths are 
the reciprocals of those altitude lengths. Therefore, crmax(^"^) < V^/ minalt(r), where minalt(r) 
means the minimum altitude. Thus, crinin(^) ^ minalt(T)/\/2. Substituting this inequality into 
H5.3[l and rearranging yields 

M - ^ ' 

where asp(r), the aspect ratio of T, equals h/ minalt(T). The aspect ratio is often used as a shape- 
quality metric; lower values mean a better shaped triangle. Thus, reversal cannot happen if the 
opposite inequality holds: 

□ 

The point of this theorem is that Type 2 reversals cannot occur for a sufficiently refined mesh 
(i.e., h sufficiently small in 15.4|l l. assuming that the mesh quality does not decay, and assuming 
that 0* is a nonsingular function. (Assuming Type 1 reversals are excluded, function 0* is never 
singular on the interior because Laplace solutions are analytic. It could be singular at the boundary 
if, for example, Q' has a corner where n had none.) 

We tested this theorem for two examples, each of which diverges a bit from the theoretical 
prediction. For the first example, we generated a uniform mesh for the rectangle [0, 2] x [0, 1] using 
Triangle and mapped all the nodes using the function /(x, y) = {x,y + ax{2 — x)). For each mesh, a 
was incremented by 1 until reversal occurred. (No Laplace solution was involved in this test case.) 
We tabulated the values of h versus a in Table[S| As predicted by the theorem, the table shows that 
as h decreases, a larger value of a is tolerated. Contrary to the theorem, however, the table shows 
that (Tmin(V/)/|| V^/ll is decreasing faster than h. In other words, reversals are avoided to a greater 
extent than predicted by the theorem. The reason for this discrepancy is that the perturbation term 
E in (|5.2|l is not well aligned with the direction that drives (V f)A toward singularity in this example. 
In particular, E affects only the y-components (since the transformation is linear for x-coordinates). 
On the other hand, transformation </) stretches the triangles substantially in the y-direction, so that 
the most effective way to perturb (V/)A toward singularity is a small change to the x-components. 

As a second test case, consider the transformation of the annulus with radii (0.5, 1) that results 
from continuous-warping continuum FEMWARP, that is, the transformation that rotates a point 
at radius p by angle a(l — r^/p^)/(l — r^), where r is the inner radius (r = 0.5 for this test). For 
each mesh, the parameter a was stepped in increments of 7r/16 until reversals were encountered. We 
tabulated values of h versus the first of a causing failure in Table As predicted by the theorem, 
decreasing h corresponded to increasing values of a, i.e., greater distortion of the domain. Again, 
these results do not initially correspond to the preceding theorem quantitively: in this case, h is 
decreasing faster than crmin(V/)/|j V^/|| . Only the last three rows of the table show that h and 
o'min(V/)/|| V^/ll are decreasing proportionally. The reason for this discrepancy is that V^/ is much 
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Table 5.1 

Second column a^^'' is the first value of a in the transformation (x, y) i— > [x,y + ax{2 — x)) of a rectangle that 
causes reversals in the mesh. The first column shows the mesh cell size (maximum edge length). The third column 
is o"jnin(V/)/||V^/|| evaluated at a vertex of a triangle that reversed. 



h 




^mi„(V/)/||V^/|| 


0.205 


10 


4.6 • 10-3 


0.108 


14 


2.5 • 10-3 


0.057 


30 


5.7- 10-* 


0.030 


51 


1.9 • 10-* 


0.015 


95 


5.6 • 10-5 



Table 5.2 

Second column a*'^" is the first value of a in the transformation (p,9) i— > {p,d(l ~ r^/p'^)/(l — r^)) (in polar 
coordinates) of an annulus that causes reversals in the mesh. The first column shows the mesh cell size (maximum 
edge length). The third column is 0"min(V/)/|| V^/|| evaluated at a vertex of a triangle that reversed. 



h 




an.in(V/)/||V^/|| 


0.202 


77r/16 


6.9- 10-3 


0.114 


tt/2 


4.6- 10-3 


0.058 


97r/16 


6.2- 10-3 


0.031 


57r/8 


2.5 • 10-3 


0.015 


ll7r/16 


2.6- 10-3 


0.008 


137r/16 


1.4- 10-3 


0.004 


TT 


7.2 • 10-* 



larger on the inner boundary than elsewhere, and the meshes in the initial rows of the table are not 
sufficiently refined to resolve the variation in the value of the derivative. 

Thus, we have seen that Type 1 reversals can be avoided by using small-step FEMWARP 
instead of FEMWARP, and Type 2 reversals can be avoided by using a sufficiently refined mesh. 
The remaining type of reverals, Type 3, are rare according to our experiments. Type 3 reversals 
are caused by the difference between the true value of the Laplace solution and the finite element 
approximation to that solution. Intuitively, this phenomenon should not be commonplace because 
the perturbation size to a mesh necessary to cause a reversal of a triangle of side- length h is 0{h), 
whereas the difference between the two mappings is 0{h?). 

Consider the following test. We generated a sequence of small-step rotations for the annulus 
in two ways. In the first case, we took steps that rotate the outer boundary by 7r/16 and leave 
the inner boundary invariant, each time computing the Laplace solution exactly analytically. This 
corresponds to iteratively applying the transformation y) = {Ax + By, —Bx + Ay) to the mesh, 
where A, B, C, D are functions of + y^ as above, and a, h, c, d are constants determined by (|5.1|l 
with r = s = 0.5, 9 = 7r/16. In the second case, we solved Laplace's equation for the above boundary 
condition using the finite element mesh that results from small-step FEMWARP. In both cases we 
tried meshes with several different values of h. The results are tabulated in TableEl As can be seen, 
discretized small-step FEMWARP outperformed continuum small-step FEMWARP. 

In other words, not only did Type 3 reversals not occur, but in fact it seems to be preferable to use 
the discretized solution for mesh warping rather than the continuum solution. The difference between 
(jf^ and 4> is the usual discretization error in finite element methods. A possible explanation for the 
improved resistance to reversals of the finite element solution is as follows. After several steps of 
small-step FEMWARP, Laplace's equation is solved on a mesh with mostly poorly-shaped elements, 
some extremely poorly-shaped. A Laplace solution minimizes the functional F{u) — Vw • Vu over 

functions u on the domain, and the finite element solution minimizes the same functional F{u) 
over the space of piecewise linear choices for u jJSl- A very poorly-shaped element is "stiffer" than 
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Table 5.3 

Continuum versus discretized small-step FEMWARP applied to a mesh of an annulus in order to test for Type 
3 reversals. The table shows that discretized small-step FEMWARP seems less prone to reversals than continuum 
small-step FEMWARP. 



h 


tail ^ 

continuum 


"discrete 


0.202 


77r/16 


7r/2 


0.114 


7r/2 


97r/16 


0.058 


97r/16 


57r/8 


0.031 


57r/8 


77r/8 


0.015 


ll7r/16 


TT 



Table 5.4 

Variable stepsize versus constant stepsize for small-step FEMWARP appied to a mesh of an annulus. Second 
and third columns are the maximum amount of rotation prior to reverals for the two methods. Fourth and fifth 
columns are the number of Cholesky factorizations required by the two methods. 



h 


„,max 




NCHOLvs 


NCHOLcs 


0.202 


1.7426 


1.7671 


13 


72 


0.114 


2.2089 


2.0862 


24 


85 


0.058 


2.6998 


2.4789 


29 


101 


0.031 


3.4852 


2.7980 


34 


114 



others in the following sense. An affine linear function u defined over a triangle that has an angle 
close to 180° will have a quite large gradient value (compared to a well-shaped triangle with the equal 
area and equal nodal values) unless the nodal values lie in a certain restricted range. Therefore, 
the extra stiffness of these elements will cause them to be deformed less than the better-shaped 
elements in the optimal solution that minimizes F. Since the poorly-shaped elements are those most 
in danger of being reversed, this is a desirable effect. 

The last topic to consider in this section is how to select a stepsize for small-step FEMWARP. 
The theory developed above indicates that as long as the step size is well below a step large enough 
to cause reversals of (p*, the step size should not matter so much. In fact, we propose the following 
simple strategy, which seems to be effective. Attempt to take a very large step (e.g., a rotation of 
size TT in the case of the annulus). If this fails (causes reversals), then halve the stepsize and try again 
until success. Update the mesh and try another such step. Note that in the process of searching for 
a correct stepsize, the coefficient matrix in FEMWARP is the same for each trial. Therefore, the 
Cholesky factors do not need to be recomputed until the mesh is updated. 

Another way to carry out small-step FEMWARP would be to take constant (small) steps on 
each iteration. We compared these two methods and found that the first was much more efficient, 
and furthermore, reversals are resisted better by the first strategy. Therefore, the repeated halving 
strategy is recommended. Table El summarizes the result for the annulus again. For the halving 
strategy, updates were pursued until the stepsize dropped below 7r/128. For the constant-step 
strategy, the stepsize was taken to be 7r/128. 

6. Mesh warping and mesh untangling. In the previous section we considered reasons why 
FEMWARP can fail and also some possible workarounds. The workarounds in the previous section, 
however, are not available in all circumstances. For example, the workaround for Type 1 reversals, 
namely, small-step FEMWARP, requires a homotopy from the old to new boundary conditions and 
also requires solution of many linear systems with distinct coefficient matrices. The workaround for 
Type 2 reversals requires refined meshes, which may not be available. 

Another workaround is to switch to a different algorithm, for example, Opt-MS, a mesh un- 
tangling method due to Freitag and Plassmann TT. Opt-MS takes as input an arbitrary tangled 
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mesh and a specification of which nodes are fixed (i.e., boundary nodes) and which are movable. 
It then attempts to untangle the mesh with a sequence of individual nodal moves based on linear 
programming. More details are provided below. 

In this section, we consider the use of Opt-MS for mesh warping. We find that the best method 
is a hybrid of FEMWARP and Opt-MS. 

To untangle the mesh, Opt-MS performs repeated sweeps over the interior nodes. For each 
interior node, it repositions the node at the coordinates that maximize the minimum signed area 
(volume) of the elements adjacent to that node (called the "local submesh"). The signed area is 
negative for a reversed element, so maximizing its minimum value is an attempt to fix all reversed 
elements. 

Let X be the location of the free vertex, that is, the current interior vertex being processed in a 
sweep. Let Xi, . . . , x„ be the positions of its adjacent vertices, and ii, . . . , t„ be the incident triangles 
(tetrahedra) that compose the local submesh. Then the function that prescribes the minimum area 
(volume) of an element in the local submesh is given by <7(x) = mini<i<„ ^^(x), where Ai is the 
area (volume) of simplex U. In 2D, the area of triangle ti can be stated as a function of the Jacobian 
of the element as follows: Ai — ^ det(xi — x, Xj — x) which is a linear function of the free vertex 
position; the same is true in 3D. Freitag and Plassmann use this fact to formulate the solution to 
max(7(x) — maxmini<i<„ ^i(x) as a linear programming problem which they solve via the simplex 
method. On each sweep, m linear programs arc solved which sequentially reposition each interior 
node in the mesh. Sweeps are performed until the mesh is untangled or a maximum number of 
sweeps has occurred. 

A shortcoming of Opt-MS, in comparison to FEMWARP, is that it is not intended to handle a 
very large boundary motion even if that motion is afhne. For example, starting from a 2D mesh, if 
the boundary vertices are all mapped according to the function [x, y) i— > (— x, —y) while the interior 
nodes are left unmoved, in many cases Opt-MS is unable to converge. FEMWARP, on the other 
hand, will clearly succeed according to Theorem 13. II 

Therefore, we propose the following algorithm: one applies first the FEMWARP mesh warping 
algorithm, and if the mesh is tangled, one uses Opt-MS to untangle the mesh output by FEMWARP 
(as opposed to the original mesh). The rationale for this algorithm is that FEMWARP is better 
able to handle the gross motions while Opt-MS handles the detailed motion better. 

To test whether this algorithm works, we applied it again to the mesh depicted in Fig. 15. II The 
boundary motion is as follows: we rotate the outer circular boundary by 9i degrees and the inner 
boundary by O2 degrees and then test three algorithms: FEMWARP alone, Opt-MS alone, and 
hybrid. (The hybrid was tested only in the case that the two algorithm individually both failed.) 
The results are tabulated in Table IHTI The table makes it clear that the hybrid often works when 
Opt-MS and FEMWARP both fail. Thus, the hybrid method is another technique for situations 
when FEMWARP alone fails and small-step FEMWARP combined with mesh refinement may be 
unavailable. The hybrid method has the disadvantage, when compared to FEMWARP, that it does 
not produce a continuous motion of interior nodes but rather only a final configuration. 

7. Three-dimensional tests. In this section we compare the robustness against reversals of 
FEMWARP, small-step FEMWARP, and the hybrid FEMWARP/Opt-MS method on examples of 
3D meshes |23 and JS]. We choose specific nonlinear boundary deformations parameterized by a 
scalar a in order to determine how much deformation each test mesh could withstand when warped 
according to each method: 

( x\ / 2 -1 \ / x\ / O.lxy \ 
\ y \ ^ \ ~2 5 y +a O.byz . 

\ z J \ Q I ) \ z j \ 0.1x2 / 

Table gives the results obtained from warping various three-dimensional meshes according to 
this boundary deformation. The data in this table is as follows. The first three columns give the 
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Table 6.1 

The row header indicates 9i (degrees), the rotation of the outer boundary, and the column header 62 is the 
rotation of the inner boundary. The table entries are as follows: 'F' means FEMWARP succeeded but not Opt-MS, 
'O' means Opt-MS succeed but not FEMWARP, 'B' means both succeeded, and 'H' means neither succeeded, but the 
hybrid succeeded, and finally '-' means none succeeded. 
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name of the mesh, the number of boundary nodes, and the number of total nodes. The next column 
•^FEMWARP is maximum value of a encountered for which FEMWARP succeeded. Parameter a 
was stepped by Aa, where Aa is indicated in the last column of the table. The table also shows 
'^ssFEMWARP' which is the last value of a for which small-step FEMWARP succeeds. The variable- 
step version of FEMWARP described in Section El was used. The minimum allowed stepsize for 
small-step FEMWARP was also taken to be Aa. 

The table indicates that small-step FEMWARP was about equal in robustness against reversals 
compared with FEMWARP except for the last mesh. In the case of "tire" , small-step FEMWARP 
was much more robust. This is probably because the reversals in the first three rows are primarily 
Type 2 reversals since the meshes are very coarse, whereas the "tire" mesh is finer and is therefore 
more likely to see Type 1 reversals according to the arguments given in the previous section. Small- 
step FEMWARP is intended to fix Type 1 reversals but is not effective against Type 2 reversals. 

We also compared Opt-MS and the hybrid FEMWARP/Opt-MS method described in the last 
section. Opt-MS performed poorly because, as mentioned earlier, it is not designed to handle large 
boundary motions. The performance of the hybrid is comparable, and in some cases superior, to 
that of small-step FEMWARP. 

Table 7.1 

Values of ct^^^ for example three-dimensional meshes. 



Mesh name 


# bdry nodes 


^ nodes 


^ max 

"FEMWARP 


max 
"ssFEMWARP 


NCHOLvs 


„ max 
"Opt-MS 


„ max 
"hybrid 


Aa 


foam5 


1048 


1337 


0.7 


1.0 


4 





1.2 


0.1 


gear 


606 


866 


3.5 


3.5 


4 


0.6 


3.5 


0.1 


hook 


790 


1190 


0.16 


0.16 


2 





0.16 


0.01 


tire 


1248 


2570 


0.15 


1.60 


13 





2.0 


0.05 



8. Application to Cardiology. We now use FEMWARP in order to study the movement of 
the beating canine heart under normal conditions. To do this, we obtained data from the Laboratory 
of Cardiac Energetics at the National Institutes of Health (NIH) [22. We were given (x, y, z, t) data 
for 192 points on the inner surface of the left and right ventricles of the beating canine heart from 
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a physiological pacing experiment. The data frames are 14.6 milliseconds apart with the first frame 
occurring 12 milliseconds before the pacing spike. 

The first step in the simulation of the ventricular movement was to generate a mesh for the 
initial position of the ventricles. In order to do this, we first noted that the 192 points we were 
given were arranged in eight slices with 24 points each. Thus, in order to generate the initial mesh, 
we decided to create a mesh for the top slice and then use FEMWARP to do the mesh warping 
necessary to create meshes for the remaining slices. Note that this uses FEMWARP in one and two 
dimensions as we describe in detail below. Once we have the meshes for all of the levels, we connect 
the triangular meshes into a tetrahedral mesh for the ventricles. The procedure to do this is also 
described below. 

We now give a more detailed description of the method we used to create the initial mesh of the 
canine ventricles. We first used Triangle to generate an initial mesh of the top slice (after projecting 
it into the x-y plane). Note that this yielded a good-quality mesh in the x-y plane with several 
additional nodes. Second, we computed the z-coordinates for the new points on the boundary of 
the top shce using ID FEMWARP. Third, we used the weight-finding portion of our FEMWARP 
algorithm to compute the weights for the appropriate 2D linear system obtained from the x- and 
y-coordinates. Fourth, we determined the z-coordinates for the mesh of the top slice by forcing the 
z-coordinates to satisfy the appropriate 3D linear system using the 2D weights. At this point, we 
had the mesh for the top slice. 

Our second task was to generate meshes for each of the remaining seven slices. This was done 
using our FEMWARP algorithm to warp the mesh for the top slice into meshes for each of the 
remaining slices. In order to accomplish this, the first step was to determine the coordinates of 
the additional boundary nodes for the mesh of the appropriate slice. This was done using ID 
FEMWARP. Then, the (x, y) coordinates of the interior nodes of that slice were determined using 
2D FEMWARP. The z-coordinates for the interior nodes were found by forcing them to satisfy the 
appropriate 3D linear system using these weights. 

The third step was to connect the triangular meshes for each of the eight slices into one tetrahe- 
dral mesh for the canine ventricles. To do this, the corresponding triangles between two slices were 
connected to form a triangular prism. After a temporary mesh of triangular prisms was created, the 
triangular prisms were subdivided into tetrahedron using the method outlined in V. 

After the initial tetrahedral mesh was created, we checked its quality using the inverse mean- 
ratio mesh-quality metric. (The mean ratio was defined in |23j and was adapted for tetrahedral 
elements in |12|.') Using this test, it was determined that the initial mesh was of poor quality. This 
was because the eight slices were equally-spaced even though the curvature of the ventricles changes 
much more rapidly near the bottom of the heart. Thus, we decided to use linear interpolation in 
conjunction with FEMWARP in the obvious way in order to add two additional slices of nodes near 
the bottom of the heart. Because the curvature of the ventricles changes more rapidly near the 
bottom of the heart, the first additional slice was placed halfway between levels 7 and 8, and the 
second additional slice was placed halfway between the first new level and level 8. The resulting 
initial mesh is shown in Figure Wa\ 

After the initial mesh was created, the heart data was used to move the 192 data points on the 
boundary of the mesh to their new positions The same process as above (i.e., FEMWARP in ID, 
FEMWARP in 2D, and using ID FEMWARP to add the two new levels of data near the bottom) 
was used in order to reposition the remaining boundary points. Once the boundary nodes were 
relocated to their positions for timestep t — 2, the 3D version of FEMWARP was used to move the 
interior nodes to their new positions for this timestep. This process was performed iteratively in 
order to study the movement of the heart at timesteps t = 3, . . . , 32. 

The simulation of the canine ventricular movement produced a series of meshes that show the 
ventricles twisting, expanding, and then contracting over the cardiac cycle. This is consistent with 
what occurs in nature. Because the dynamic range of the motion is small, it cannot be detected in 
single figures (separate from an animation). 
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Fig. 8.1. Initial canine ventricular mesh. 

During each timestep, the inverse mean ratios of the tetrahedra in the mesh were computed. The 
inverse mean ratio computations showed that the heart remained untangled throughout the entire 
simulation. This is not very surprising since the heart mesh is composed of elliptical rings that seem 
to undergo less movement on each timestep than the circular rings in the test cases. However, the 
motion of the heart is anisotropic which makes it difficult to predict in advance how it will tolerate 
deformations. Interestingly enough, the values of the minimum and average inverse mean ratios were 
relatively constant across all timesteps. Only the value of the maximum inverse mean ratio changed 
a significant amount. However, it only increased by as much as four percent of its initial value and 
decreased by as much as seven percent of that value. Thus the inverse mean ratio computations 
indicate that the heart meshes are of sufficiently good quality for use with a numerical PDE solver 
that requires moving meshes. 

In order to further test our mesh warping algorithm, the motion of the ventricles was exaggerated 
by a factor of 3. In this case, the motion was large enough to detect in separate figures and is shown 
in Figure 18.21 We note that the FEMWARP algorithm also performed successfully in this case, 
which is encouraging given the much larger deformations. Small-step FEMWARP was not needed 
for this problem. 

9. Conclusions. We studied an algorithm called FEMWARP for warping triangular and tetra- 
hedral meshes. The first step in the algorithm is to determine a set of local weights for each interior 
node using finite element methods. Second, a user-supplied deformation is applied to the boundary 
notes. The third and final step is to solve a system of linear equations based upon the weights and 
the new positions of the boundary nodes to determine the final positions of the interior nodes. 

FEMWARP falls into a more general class of methods that we call linear weighted laplacian 
smoothing (LWLS). LWLS methods have several advantages. First, if a continuous boundary defor- 
mation is given, then methods within the framework are valid for computing the resulting trajectory 
that specifies the movement of the interior nodes. In addition, these trajectories will be continuous, 
which is vital for applications where continuity of motion is required. Second, sparse matrix algo- 
rithms may be used to solve the linear system which determines the final positions of the interior 
nodes. Third, the methods are exact if the boundary deformation is afhne. 

We then analyzed the case that FEMWARP fails, i.e., produces element reversals. Some 
workarounds include: taking smaller steps, using a finer mesh, and finally, a hybrid of FEMWARP 
and Opt-MS. 

We tested the robustness of FEMWARP, small-step FEMWARP, and hybrid FEMWARP/Opt- 
MS on 2D annulus test cases and 3D general unstructured meshes. The latter two algorithms 
generally outperform plain FEMWARP, sometimes significantly. We also used FEMWARP to study 
the motion of the canine ventricles under normal conditions. 
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